,  AD-A130  297 

A  CHARACTERIZATION  OF  ARCTIC  UNDERSEA 
PRINCETON  UNIV  NJ  DEPT  OF  STATISTICS 

NOISE(U) 

J  G  VEITCH  ET  RL. 

1/1 

N 

JUN  S3  TR-12  N80014-81-K-0146 

UNCLASSIFIED 

F/G  14/2 

NL 

IES01UTI 

fE.AU  OF  ! 


1C  FILE  COPY 


REPORT  NUMBER  12 


ADA  I  3<>39;; 


A  CHARACTERIZATION  OF  ARCTIC 
UNDERSEA  NOISE 


J.G.  VEITCH  and  A.R.  WILKS 


Department  of  Statistics 
Princeton  University 
Princeton,  New  Jersey  08544 


JUNE  1983 


Prepared  for 

OFFICE  OF  NAVAL  RESEARCH  (Code  41  ISP) 

Statistics  and  Probability  Branch 

Arlington,  Virginia  22217 

under  Contract  N00014-81-K-0146 

SRO(103)  Program  in  Non-Gaussian  Signal  Processing 

S.C.  Schwartz,  Principal  Investigator 


Approved  for  public  release;  distribution  unlimited 

83  07  14  002 


^ELEC-nfj 

JUL  1  4  1983 


«ECU«ITY  CLASSIFICATION  OF  THIS  PAGE  fPTim  Om  (ni»n« 


REPORT  DOCUMENTATION  PAGE 


I.  REPORT  NUMBER 


4.  TITLE  ran*  Suktlllm) 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


S.  RECIPIENT'S  CATALOG  NUMBER 


s.  type  of  report  a  period  covered 

Technical  Report 

Sept.  1982-May  1983 

s.  performing  org.  REPORT  number 


a.  contract  or  grant  numbers; 


N00014-8 l-K-014  6 


10.  PROGRAM  ELEMENT.  PROJECT.  Task 
AREA  •  PORK  UNIT  NUMBERS 


A  CHARACTERIZATION  OF  ARCTIC  UNDERSEA 
NOISE 


7.  AUTHORS  FI 

James  G.  Veitch  and  Allan  R.  Wilks 


t.  performing  organization  name  ano  address 

Department  of  Statistics 
Princeton  University 
Princeton,  New  Jersey  08544 


•  I.  CONTROLLING  office  NAME  and  aooress 

Office  of  Naval  Research  (Code  411  SP) 
Department  of  the  Navy 
Arlington,  Virginia  22217 


i«  MONITORING  AGENCY  name  A  ADDRESS*’!/  Clllaraai  Bam  Cantralllat  Oltlea)  IS.  SECURITY  CLASS.  (•!  till*  tapan) 

Unclassified 


NR  SRO-103 


12.  REPORT  OATS 

June  1983 


IS.  NUMBER  OF  PAGES 

42 


1ft.  DISTRIBUTION  STATEMENT  (of  U%lo  Roport) 


Approved  for  Public  Release;  Distribution  Unlimited 


17.  DISTRIBUTION  STATEMENT  (at  Uta  akatrael  MirM  In  Black  it.  It  ttflaranl  Bam  ftapprlj 


<S.  KEY  POROS  (Caniinaa  an  r araraa  alta  It  naaaaaarr  nrP  Itaniltr  A y  Hack  w— Air, 


arctic  noise 

noise  analysis  and  models 

arctic  noise  spectrum  analysis 


20.  ABSTRACT  rCpilliiw  m  nwh  cl  am  It  nmaaamaky  m4  Ikantlty  Ar  Hack  mamkat) 

Arctic  undersea  noise  is  investigated  and  partially  character- 
zed  by  a  probability  model  suggested  by  exploration  of  a  sample  of 
uch  noise.  The  large  size  of  the  dataset  makes  this  possible, 
he  model  is  that  of  a  mixed  spectrum,  with  stationary  Gaussian 
oise  intermixed  with  random  sinusoids,  and  occasional  high  inten¬ 
sity  impulsive  noise  bursts,  which  are  very  short-lived.  Due  to  the 
complex  behavior  of  the  random  sinusoids,  the  noise  cannot  be  con¬ 
sidered  as  stationary  in  the  time  frame  (10  minutes)  investigated. 


I  JAN  71 


EDITION  OF  I  NOV  ••  IS  OBSOLETE 

S/M  01 02*  LF  4!  44401 


(con't.  ) 

SECURITY  CLASSIFICATION  OF  TMIS  page  (PR  IN  Data 


_  I 

SCCUNtTV  CLASSIFICATION  OF  THIS  PAOC  fWhwi  Data  Cnlara« 


(Abstract  con't.) 

The  validity  of  the  model  is  checked  by  testing  model  predictions 
against  the  actual  data.  In  particular,  this  model  makes  some 
predictions  about  the  type  of  spectrum  observed,  and  these 
predictions  are  tested  by  a  non-linear  filtering  of  the  period- 
ogram.  A  running  median  filter  is  used  and  is  shown  to  asymptot¬ 
ically  estimate  the  spectrum  of  the  Gaussian  noise.  This  filter 
seems  to  work  fairly  well  on  the  data.. 


sicuNirr  classification  of  tnis  faowwmn  »«*•  *»*•»•« 


:NSPi.CTftl 

V  3  V 


A  Characterization  of  Arctic  Undersea  Noise 

by 

James  G.  Veitch  and  Allan  R.  Wilks 
Department  of  Statistics 
Princeton  University 
Princeton.  New  Jersey,  06544 

Abstract 

Arctic  undersea  noise  is  investigated  and  partially  characterized  by  a  pro¬ 
bability  model  suggested  by  exploration  of  a  sample  of  such  noise.  The  large 
size  of  the  dataset  makes  this  possible.  The  model  is  that  of  a  mixed  spectrum, 
with  stationary  Gaussian  noise  intermixed  with  random  sinusoids,  and  occasional 
high  intensity  impulsive  noise  bursts,  which  are  very  short-lived.  Due  to  the 
complex  behavior  of  the  random  sinusoids,  the  noise  cannot  be  considered  as 
stationary  in  the  time  frame  (10  minutes)  investigated.  The  validity  of  the 
model  is  checked  by  testing  model  predictions  against  the  actual  data.  In  par¬ 
ticular,  this  model  makes  some  predictions  about  the  type  of  spectrum 
observed,  and  these  predictions  are  tested  by  a  non-linear  filtering  of  the 
periodogram.  A  running  median  filter  is  used  and  is  shown  to  asymptotically 
estimate  the  spectrum  of  the  Gaussian  noise.  This  filter  seems  to  work  fairly 
well  on  the  data. 
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1.  Description  of  Data  and  Goals  of  Analysis 

This  report  describes  the  analysis  of  a  particular  dataset.  In  this  section  we 
describe  the  dataset,  and  the  goals  of  the  analysis.  In  section  2  we  look  at  some 
previous  work  done  by  others  in  analyzing  this  and  other  similar  datasets.  In 
the  remaining  sections  we  describe  the  steps  in  our  data  analysis,  along  with 
some  theoretical  justification  of  severed  novel  techniques. 

The  dataset  is  a  large  univariate  time  series.  It  is  a  sample  of  underwater 
noise  observed  on  23-24  April  1980  from  a  pack  ice  camp  in  the  Arctic  Ocean, 
located  at  86  degrees  north  latitude,  25  degrees  west  longtitude.  A  broadband, 
omnidirectional  hydrophone  was  suspended  to  a  depth  of  91  meters  in  4000 
meter  deep  water,  and  the  underwater  acoustic  noise  was  recorded  on  an  analog 
device.  Later,  this  noise  was  filtered  by  a  lowpass  filter  with  cutoff  around  2500 
Hz,  and  then  digitized  at  10  kHz,  ie,  the  analog  signal  was  sampled  at  the  rate  of 
10,000  times  per  second,  and  each  of  these  sampled  values  was  rounded  to  the 
nearest  of  a  finite  set  of  numbers.  Our  dataset  is  this  digitized  record  during 
the  period  approximately  11:30-11:40  pm  on  the  23rd.  Ten  minutes  at  10  kHz 
yields  six  million  numbers:  we  actually  have  6,150,144  numbers  which  is  slightly 
longer  than  ten  minutes.  This  conveniently  works  out  to  6006  groups  of  1024 
numbers  each. 

Such  a  large,  single  time  series  presents  both  problems  and  opportunities; 
this  is  perhaps  a  general  comment  relevant  to  the  analysis  of  large  datasets. 
The  main  problem  is  that  the  statistics  we  usually  compute  such  as  moments, 
autocovariances  and  Fourier  transforms,  are  extremely  inconvenient  to  com¬ 
pute  for  large  amounts  of  data,  and  may  even  make  little  sense.  Just  handling 
the  data  presents  logistic  problems.  On  the  other  hand,  we  can  do  things  which 
are  impossible  or  impractical  for  smaller  datasets.  Many  such  techniques  will  be 
described  in  the  following  pages  and  perhaps  what  makes  these  techniques 
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attractive  is  the  general  idea  that  by  blocking  the  data  into  contiguous  chunks 
we  caui  imagine  that  we  have  an  ensemble  of  time  series.  This  contrasts  with  the 
usual  case  in  time  series  analysis  in  which  we  only  have  one  time  series  from  a 
conceptual  ensemble  of  possible  time  series. 

The  particular  dataset  analyzed  here  appears  to  exhibit  several  different 
sorts  of  non-stationarity  and  non-Gaussianity.  Our  main  goal  in  analyzing  this 
so-called  'ice  noise’  is  to  try  to  characterize  it.  We  wish  to  find  a  model  which 
adequately  describes  the  noise.  This  characterization  could  be  used  to  enhance 
knowledge  about  underwater  acoustics  near  ice,  and  might  also  allow  more 
efficient  detection  of  signals  in  the  presence  of  such  noise. 

We  have  taken  a  data-oriented  approach  in  our  analysis.  We  begin  by  pic¬ 
turing  the  data  in  certain  simple  ways  and  then  letting  these  pictures  guide  our 
further  anedysis.  The  end  result  is  an  approximate  theoretical  model,  suggested 
to  us  purely  by  consideration  of  the  data  itself.  It  is  to  be  expected  that 
detailed  knowledge  of  underwater  acoustics  under  ice  would  lead  to  at  least 
some  enhancement  of  our  model.  This  point  will  be  taken  up  later  in  greater 
detail. 

2.  Previous  Work 

This  dataset  has  been  looked  at  several  times  before.  In  particular.  Roger 
Dwyer  (1981,  1982),  has  done  considerable  work  in  which  he  investigates  running 
moments  in  the  time  domain  and  does  analysis  in  the  frequency  domain.  He 
concludes  that  the  noise  is  both  non-Gaussian  and  nonstationary.  In  particular 
he  sees  very  large  sample  kurtoses,  which  are  attributed  to  impulsive  noise,  and 
notes  the  presence  of  pronounced  tonals  in  the  spectrum  estimates  which  are 
hypothesized  to  be  squeaks  due  to  ice  dynamics  Similar  results  have  been 
reported  by  Milne  and  Ganton  (1964)  in  other  data  sets 
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we  confirm  the  presence  of  such  impulsive  noise  and  such  tonals,  but 
believe  that  the  remainder  of  the  under  ice  noise  is  essentially  Gaussian  in 
nature.  This  is  explored  in  Section  4. 

3.  Data  Blocking  and  Some  Simple  Traces 

Many  of  the  commonly  used  techniques  in  time  series  rest  on  the  twin 
assumptions  of  stationarity  and  Gaussianity.  The  largeness  of  our  dataset  gives 
us  the  opportunity  to  examine  these  assumptions  for  validity.  Since  computing 
some  of  the  traditional  statistics  such  as  moments,  quantiles  and  Fourier 
transforms  is  impractical  for  six  million  numbers,  we  are  naturally  led  to  the 
idea  of  blocking  the  data.  To  do  this  we  group  consecutive  observations  into 
groups  of  size  N,  say,  where,  if  convenient,  N  may  be  taken  as  a  divisor  of  T,  the 
total  number  of  data  points.  By  computing,  for  example,  moments  for  each  of 
these  R  =  77  TV  records,  we  may  ascertain  whether  stationarity  and  Gaussianity 
are  reasonable  assumptions. 

3. 1.  Choosing  a  Block  Size 

From  physical  considerations,  there  is  no  natural  record  size  N,  for  our  ten 
minutes  of  acoustic  noise,  and  we  are  thus  free  to  choose  it  as  we  please.  In 
most  of  our  analysis  we  have  used  N  =  1024  for  the  following  reasons: 

(i)  Dwyer  has  used  this  value  in  his  analyses,  so  our  work  and  his  may  be  com¬ 
pared  at  appropriate  points. 

(ii)  A  power  of  two  is  convenient  for  the  computation  of  discrete  Fourier 
transforms  via  the  FFT  algorithm. 

(iii)  ,V  should  be  greater  than  about  100  in  order  to  achieve  satisfactory  stan¬ 
dard  errors  for  record  statistics,  and  less  than  about  10,000  in  order  that 
the  computations  do  not  get  too  cumbersome.  Our  chosen  N  is  roughly  the 
geometric  mean  of  these  two  extremes. 


(iv)  N  should  be  sufficiently  small  that  within  a  single  record  the  data  is  roughly 
stationary.  1024  seems  empirically  to  satisfy  this  condition,  though  more 


rigorous  testing  would  be  possible. 

(v)  N  should  be  at  least  512  since  this  blocksize  corresponds  to  .05  seconds  of 
data  so  that  the  frequency  of  blocks  is  20  Hz,  and  signals  of  10  Hz  or  less 
would  appear  as  a  trend  from  block  to  block.  Such  frequencies  are  unlikely 
to  be  detected  with  the  recording  equipment  used,  though  with  a  smaller 
blocksize,  there  might  start  to  be  block  to  block  trends. 

Thus,  our  first  step  is  to  break  the  data  into  R  =  6006  records  each  of  size 

N  =  1024  numbers. 


3.2.  Running  Means  and  Variances 

Our  first  plots  replicate  those  in  Dwyer  (1901),  page  0  Denote  the  jth 
record  by  xi  0,  .  .  ,  for  j  -  1 ...  ,R,  with  N  -  1024  and  R  -  6006.  Thus  the 

entire  time  series  is  just  the  numbers 


•  >  *  i.io23'  •  .*aoc6.n . *ecca.;oz3j 

We  compute  for  the  ;th  record  the  mean 


IJS3 


m,  = 


1024  l 


i*0 


‘■j.i 


and  the  sample  central  moments 


k  =  2,3,4. 


Each  of  these  sample  moments  is  computed  for  records  which  do  not  con¬ 
sist  of  independent  observations,  but  in  fact  are  highly  correlated.  Thus  the 
sample  moments  will  be  biased  estimates  of  the  true  moments.  Notice  that  we 
are  tacitly  assuming  that  the  data  is  relatively  stationary  within  records,  so  that 
it  makes  sense  to  speak  of  the  theoretical  moments  of  a  record.  As  an  example 
of  this  bias,  consider  the  sample  variance  in  the  presence  of  positive  serial 
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correlation.  In  this  case  sample  variance  would  estimate  the  variance  with  a 
negative  bias,  i.e.  it  would  underestimate  the  true  variance.  We  choose  to  ignore 
the  effects  of  bias  in  what  follows  because  the  sample  size  is  relatively  large 
(1024)  and  because  the  records  tend  to  have  sample  autocorrelation  functions 
which  drop  off  fairly  rapidly 

The  actual  center  and  variability  of  the  raw  data  are  artifacts  of  the  digiti¬ 
zation  method  used  on  the  original  analog  signal.  Thus,  we  denote  by  m  the 
sample  average  of  the  entire  time  series,  and  by  the  sample  variance  of  the 
time  series,  and  plot  in  Exhibits  3.1  and  3.2  the  values  of  (m.j  -  m)/  (m^)1/2 
and  m/aV  m(2),  as  functions  of  j , 

It  is  clear  from  Exhibit  3.1  that  the  data  seem  trendless,  and  that  except 
for  records  2220,  2248  and  3236,  we  are  just  seeing  sample  fluctuations  in  the 
sample  average.  This  trendlessness  is  to  be  expected,  since  the  hydrophone  was 
measuring  longtitudinal  pressure  waves  which  must  clearly  average  out  to  zero 
in  any  length  of  time  relatively  long  compared  to  the  period  of  the  longest  wave 
present  (See  reason  (v),  above,  for  the  use  of  a  blocksize  of  1024.) 

The  sample  variance  trace,  Exhibit  3.2  clearly  demonstrates  that  over  the 
time  period  of  the  data  there  is  marked  non-stationarity.  In  fact  it  would  appear 
that  even  in  the  short  run  there  is  no  stationarity,  with  the  possible  exception  of 
a  few  sections  of  the  data  where  the  sample  variance  trace  is  relatively  flat. 
Sample  variance  is  a  measure  of  volume  or  signal  energy,  and  thus  it  appears 
that  the  volume  of  the  noise  is  fluctuating  In  later  analysis  we  will  show  that 
much  of  this  variability  in  the  variance  can  be  attributed  to  the  momentary 
presence  of  strong  sinusoidal  components  superimposed  upon  relatively  stable 
noise,  and  that  the  variance  fluctuation  is,  in  fact,  not  being  caused  by  changes 
in  the  volume  of  this  'background'  noise. 
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3.3.  Running  Skews  and  Kurtoses 

Exhibits  3.3  and  3.4  are  plots  of  the  traces  of  the  sample  skew 

mfs) 

1  (mpf/2 ' 

and  the  sample  kurtosis 

K.=  ^ 

for  each  of  the  6006  records, 

The  skew  and  kurtosis  for  any  Gaussian  distribution  are  0  and  3.  Exhibits 
3.3  and  3.4  show  skews  and  kurtoses  which  are  usually  reasonable  for  a  Gaussian 
distribution,  with  occasional  dramatic  deviations  from  the  baselines  of  0  and  3. 
In  looking  at  these  two  plots  together,  it  is  striking  that  the  large  deviations  in 
skew  and  kurtosis  appear  to  occur  together  much  of  the  time. 

3.4.  Explanation  for  Large  Skews  and  Kurtoses 

To  illustrate  this  point  the  absolute  skews  and  the  kurtoses  were  each 
sorted,  and  the  50  largest  selected  from  each  list.  Those  records  which  are 
common  to  the  two  lists  are  sorted  by  kurtosis  in  Table  1  below.  Nearly  half  of 
the  spikes  in  the  skew  and  kurtosis  traces  occur  for  the  same  records  among 
the  largest  50  of  each.  In  addition,  the  largest  12  absolute  skews  and  the  largest 
1 1  kurtoses  appear  on  this  list.  Thus  large  absolute  skew  and  large  kurtosis  are 
highly  correlated.  Exhibit  3.5  is  a  time  series  plot  of  the  data  for  record  1362, 
which  has  the  second  highest  kurtosis  and  the  largest  absolute  skew.  The  rea¬ 
son  for  the  large  moments  is  immediately  clear,  viz.  the  short  burst  of  very 
large  observations  around  observation  number  695.  The  duration  of  the  burst  is 
about  10  observations  which  in  real  time  is  about  1  millisec  ond  This  behavior  is, 
in  fact,  characteristic  of  records  with  larg  lew  iQ  large  kurtosis.  and 
accounts  for  large  kurtoses  at  least  as  small  as  6,  as  illustrated  by  Exhibit  3.6,  a 


time  series  plot  of  record  2238.  This  record  has  the  50th  largest  kurtosis  (5  73) 
and  a  relatively  small  skew  (.39). 


Table  1 

Largest  Skews  and  Kurtoses 

Record  Kurtosis  Skew 


2066 

29.51  -1.57 

1362 

28.17  -3.12 

4947 

27.91  -2.55 

1487 

24.15  0.51 

2220 

22.97  -2.72 

2248 

18.49  -2.92 

41 

17.67  -0.68 

3236 

17.07  -2.76 

5776 

14.95  0.62 

2236 

14.54  -1.39 

2177 

13.41  -1.05 

1281 

12.42  -0.70 

2041 

11.25  -0.89 

2246 

10.24  1.02 

4902 

10.00  -1.04 

5546 

9.99  -0.67 

5557 

9.55  -0.84 

5705 

9.25  -0.64 

4808 

8.05  -0.52 

2042 

7.43  0.52 

1938 

7.10  -0.55 

2724 

6.59  -0.63 

4836 

6.12  -0.66 

Kurtosis  is  classically  used  to  measure  "fatness"  in  the  tails  of  a  distribution 
and  these  bursts  appears  to  be  the  major  factor  contributing  to  the  high  Kn  If. 
for  sample  1362,  we  excise  the  points  from  689  to  699  and  recalculate  kurtosis, 
Kl3S2  drops  from  28.17  to  3.24.  Skew  is  also  significantly  affected  by  these 
bursts;  for  sample  1362,  the  same  excision  causes  S i362  to  go  from  -3.12  to 
0.039. 

As  stated  previously,  kurtosis  measures  "fatness"  in  the  tails  of  the  distri¬ 
bution.  However,  a  large  kurtosis  does  not  tell  us  how  far  from  the  center  of  the 
distribution  we  must  go  before  we  see  "fat”  tails.  We  wished  to  check  that  the 
observed  high  sample  kurtoses  were  being  caused  by  bursts  (i.e.  the  tails  are 


"fat”  only  far  from  the  center).  For  we  expect  that  the  body  of  the  distribution 
of  the  numbers  from  a  record  will  remain  unchanged  by  a  burst,  but  that  the  few 
largest  and  smallest  values  in  the  sample  will  be  much  larger  (or  smaller, 
respectively).  We  devised  a  check  by  computing  the  statistic: 

Fn  =  Yn/Hn 

where  Yn  is  the  difference  between  the  second  highest  and  second  lowest  obser¬ 
vations  and  Hn  is  the  interquartile  range. 

This  statistic  is  an  analog  of  the  sample  kurtosis,  in  that  it  measures  how 
long  the  tails  are  ( Yn )  relative  to  the  spread  of  the  distribution  ( Hn ).  Sample 
kurtosis  achieves  this  effect  by  dividing  the  average  fourth  power  of  deviations 
from  the  mean  by  the  average  second  power.  The  latter  measures  spread,  while 
the  former  causes  greater  magnification  of  those  observations  which  deviate 
further  from  the  mean.  The  problem  with  using  kurtosis  as  a  measure  of  the 
length  of  the  tails  of  a  distribution  is  that  kurtosis  can  also  be  large  for  a  distri¬ 
bution  which  is  peaked,  i.e.  a  distribution  whose  tails  are  fat  relative  to  the 
center,  but  whose  extreme  tails  are  not  particularly  large.  Fn  will  only  be  large 
in  the  presence  of  large  extreme  tails  in  the  distribution  and  can  therefore  be 
used  to  check  our  belief  that  the  observed  large  sample  kurtoses  for  the  ice 
noise  are  due  to  very  high  amplitude  bursts  and  not  merely  to  a  peaked  distri¬ 
bution. 

For  a  comparison  of  the  two  measures  see  Exhibit  3.7,  which  illustrates  the 
trace  of  Fn ,  limited  to  between  2  and  7,  and  the  trace  of  limited  to  between  2 
and  5.  The  values  were  so  limited  in  order  that  the  behaviour  of  the  majority  of 
the  values  was  not  swamped  by  the  few  large  values  in  each  trace.  It  is  clear 
that  the  two  traces  are  virtually  the  same  in  character.  The  correlation  between 
the  logs  of  the  measures  is  .93  (The  reason  for  taking  logs  is  to  reduce  the 
effect  the  few  large  values  have  on  the  correlation  coefficient.)  The  group  of 
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records  with  the  127  highest  kurtoses  and  the  group  of  records  with  the  127 
highest  Fn's  share  117  members  in  common!  Kurtosis  and  Fn  measure  essen¬ 
tially  the  same  thing,  and  therefore  kurtosis  is  large  precisely  because  of  the 
stray  values  represented  by  bursts. 

A  fair  question  to  ask  is  the  following:  If  we  wish  to  automatically  detect 
bursts  in  this  noise,  how  large  should  the  kurtosis  of  a  record  be  before  we  flag 
such  a  burst?  We  settled  on  4,  after  considering  a  portion  of  the  data  analyzed 
in  section  4.1.1.  In  this  portion  of  the  data  all  but  three  kurtoses  are  under  3.9; 
the  remaining  three  are  all  much  larger.  Theoretic al  considerations  also  make 
this  a  generous  upper  bound  -  see  section  4.1.1  for  a  discussion  of  simulation 
results. 

It  would  appear,  then,  that  all  the  large  kurtosis  records  have  a  common 
cause:  a  transient  large  amplitude  signed  which  disappears  very  rapidly  com¬ 
pared  to  the  time  scale  of  each  record.  The  erratic  behavior  of  this  transient 
seems  to  account  for  the  large  skew  values  often  observed  along  with  large  kur¬ 
toses.  Because  Kn  demarcates  these  records,  we  have  a  simple  method  for  iden¬ 
tifying  the  records  containing  bursts. 

3.5.  Other  Behavior  in  the  Skews  and  Kurtoses 

Since  these  bursts  are  so  large,  they  obscure  other  behavior  in  the  moment 
traces  for  skew  and  kurtosis  in  Exhibits  3.3  and  3.4.  In  order  to  see  this  other 
behavior,  we  repiot  in  Exhibit  3.8  the  skew  and  kurtosis  traces  in  the  regions  of 
the  0  and  3  baselines.  The  variance  trace  is  also  replotted  in  this  Exhibit.  This 
Exhibit  contain  most  of  the  skews  and  kurtoses;  in  fact  99.6%  of  the  skews  and 
98.8%  of  the  kurtoses.  In  the  skew  plot,  there  is  a  region  of  relatively  large  nega¬ 
tive  skews  from  about  300  to  1400  and  two  segments  where  the  skews  seem  to  be 
compressed  toward  the  origin.  1600-1900  and  3300-4300  In  the  kurtosis  plot 
there  seem  to  be  three  regions  where  the  kurtosis  drops  below  3:  1000-1400, 
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1600-1900  and  3300-4300.  These  regions  also  seem  to  form  natural  boundaries 
on  the  variance  plot.  The  coincidences  of  these  regions  bears  investigation;  we 
will  suggest  a  model  in  the  next  section  which  appears  to  account  for  these 


effects.  To  this  end,  we  divide  and  name  the  data  as  follows; 

Region  Records  Variance  Skew  Kurtosis 


I 


1 

2 

3 

4 

5 

6 

7 

8 


1-300  normal  normal 
301-1000  high  stretched 

1001-1400  high  stretched 

1401-1600  normal  normal 


1601-1900  high  squeezed 
1901-3300  normal  normal 
3301-4300  high  squeezed 
4301-6006  normal  normal 


normal 

normal 

low 

normal 

low 

normal 

low 

normal 


4.  The  Basic  Model. 


On  inspection  of  the  variance  trace,  it  is  clear  that  the  volume  of  the  noise 
appears  relatively  constant  over  fairly  long  stretches  of  time.  One  can  see  this 
^  on  casual  inspection  of  the  variance  trace  in  Exhibit  3  2,  for  record  numbers 

1900-3300  or  4300-6006  (regions  6  and  8).  In  the  regions  where  the  variance  is 
fluctuating  much  more  wildly,  many  of  the  log  periodograms  of  single  records 
I  have  very  marked  peaks  at  a  few  frequencies.  Further,  for  these  log  periodo¬ 

grams,  on  ignoring  peaks  "by  eye",  the  underlying  periodogram  looks  very  much 
the  same  from  record  to  record  This  suggests  that  the  ice  noise  might  be  well 
•  modeled  by  a  continuous  spectrum,  stationary  series  with  added  discrete 

sinusoids  of  high  power  whose  frequencies  change  from  record  to  record.  An 
added  complexity  would  be  to  allow  the  continuous  portion  to  vary  slowly 
i  through  time  We  shall  call  such  discrete  sinusoids  tones,  or  tonais  Mathemati¬ 

cally,  let  F„(cj)  be  the  observed  periodogram  of  record  n;  u  ranges  through  the 
Fourier  frequencies.  Then  the  former  hypothesis  is  essentially  one  of  a  mixed 
,  spectrum: 


Pn{u)  =  h(u)  +  £  A^Fiiu-ipin)  +  sn(u) 


(1) 


The  first  term  on  the  right  is  the  continuous  portion  represented  by  a  spectral 
density  function  h(cS).  The  sum  represents  discrete  tones  at  frequencies 
V\n.  Via  f°r  some  small  integer  /.  and  presents  the  effect  of  computing  a 
periodogram  as  a  sum  Fejer  kernels  Fu...,Fi  weighted  by  the  amplitudes 
•dm.  .Ajn  of  the  tones  If  we  allow  variation  of  the  continuous  part,  the  model 
becomes: 

P n(^)  =  ^n(^)  +  ^  AnFi(w— Pin)  £n(y)  (2) 

i=l 

where  the  may  be  slowly  varying  with  n  in  a  complex  fashion.  The  tones  are 
allowed  to  be  quite  arbitrary  from  record  to  record;  the  only  constraint  is  that 
there  be  no  more  than  I  of  them  at  any  one  time.  If  we  assume  in  both 
hypotheses  that  the  part  of  the  process  with  spectrum  h(u)  (or  ^(u),  in  Model 
(2))  is  stationary  Gaussian,  then  the  cn(«)  are  independent,  and  exponentially 
distributed  with  mean  0  and  variance  h(6).  (Note  that  the  cn(v)  are  exponential 
because  Pn  is  the  raw  periodogram;  it  is  the  sum  of  the  squares  of  two  indepen¬ 
dent  normals  -  the  real  and  imaginary  parts  of  the  discrete  Fourier  transform) 
We  shall  call  the  component  of  the  noise  with  spectral  density  h  the  background 
noise 

We  shall  further  assume  (as  is  suggested  by  the  log  periodograms  that  we 
looked  at)  that  h  (or  hn )  is  fairly  smooth.  In  particular,  we  assume  that  any 
peaks  in  h  are  much  wider  than  the  width  of  the  peaks  of  tones  appearing  in  the 
periodogram. 

Both  of  these  models  have  some  flaws  that  are  apparent  on  casual  inspec¬ 
tion.  First,  the  high  frequency  part  of  all  the  raw  periodograms  are  very  flat 
(e.g  Exhibit  4.1  pictures  a  periodogram  from  record  4601)  -  the  error  distribu¬ 
tion  e*(u)  is  obviously  only  poorly  modeled  here.  This  high  frequency  part  of  the 
spectrum  has  such  low  power  that  modifying  any  single  observation  can  change 
the  high  frequency  portion  of  the  periodogram  by  several  orders  of  magnitude. 


We  should  recall  that  the  data  was  filtered  by  a  low-pass  filter  with  a  cut-off  at 
2500  Hz,  so  it  is  not  surprising  the  models  fail  in  this  region.  We  shall  thus 
ignore  the  frequencies  above  2500  Hz.  Second.  Model  (!)  appears  flawed  in  that 
even  in  the  absence  of  tones  h  seems  to  vary  systematically  over  time  --  look  at 
the  variance  trace  for  records  4430-4450.  That  there  are  no  discrete  frequen¬ 
cies  for  these  records  will  be  further  explored  in  the  section  on  background 
noise,  Section  4.1.1.  It  is  clear  that  any  such  model  as  (1)  can  be  at  best  only 
approximately  true. 

There  is  a  third  flaw  with  Model  (1)  In  some  of  the  records  very  large  tran¬ 
sient  pressure  changes  are  obvious,  as  discussed  in  Section  3.4.  These  transient 
bursts  are  on  the  order  of  1/100  of  a  second  or  less.  They  seem  to  occur  in  only 
a  small  fraction  of  the  records.  We  could  modify  Model  (l)  or  (2)  to  accommo¬ 
date  these  bursts  by  adding  a  "burst"  to  the  process  at  random  intervals,  but  as 
records  with  bursts  can  be  identified  as  discussed  in  section  3.4  we  shall  assume 
that  conclusions  about  Model  (l)  or  (2)  apply  to  these  records  (except  for 
bursts),  and  we  shall  base  our  study  of  Model  (l)  on  records  without  these 
bursts.  We  do  not  attempt  to  model  bursts  in  this  analysis. 

Our  aim  is  to  check  whether  Model  (1)  is  a  reasonable  description  of  the 
data,  and  this  will  be  done  in  the  following  sections. 

4.1.  Confirmation  of  Model  (l) 

4. 1. 1.  Gaussianity  in  the  Background  Noise. 

We  can  reasonably  assume  there  are  no  tones  in  areas  where  the  variance 
trace  is  flat.  We  thus  wish  to  check  whether  the  noise  in  these  regions  is  station¬ 
ary  Gaussian;  this  would  be  predicted  by  Model  (1)  In  order  to  substantiate  this 
assumption  we  considered  the  sequence  of  records  4601-4901.  A  log  periodo- 
gram  of  one  of  these  (4601)  is  given  in  Exhibit  4.1.  We  ignore  three  of  these 
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records,  viz  4692,  4710  and  4808,  where  large  kurtoses  of  6.38,  11  36  and  8.05 
flagged  bursts  in  the  noise.  Otherwise,  all  the  moment  traces  are  very  flat,  with 
the  variance  trace  ranging  between  378  and  .764.  The  observed  kurtoses 
remain  close  to  three  with  an  average  of  3.02  and  a  standard  deviation  of  .233. 
In  a  simulation  of  10,000  sample  kurtoses  from  records  of  1024  simulated 
independent  Normals  we  obtained  an  average  kurtosis  of  2.995,  with  a  standard 
deviation  of  .152;  this  is  not  strictly  comparable  since  we  have  dependence  in 
the  ice  noise,  but  it  suggests  that  the  ice  data  kurtoses  are  in  fact  what  one 
might  expect  under  an  assumption  of  Gaussianity.  We  get  very  similar  results 
for  skew;  the  ice  data  skews  for  these  records  average  .017  with  a  standard  devi¬ 
ation  of  .125,  to  be  compared  with  an  average  of  -.00007,  and  a  standard  devia¬ 
tion  of  .0785  for  the  simulated  data.  This  is  encouraging  evidence  that  records 
4601  to  4901  may  be  Gaussian  noise. 

We  substantiate  this  further  by  considering  the  actual  distribution  of  the 
data  for  part  of  this  region.  However,  the  difficulty  of  handling  so  many 
numbers  at  once  means  we  shall  do  this  for  only  a  small  proportion  of  these 
numbers  The  technique  we  use  is  to  get  a  large  independent  sample  by  taking 
every  64th  number  in  this  region.  This  is  reasonable  since  autocorrelations  tend 
to  fluctuate  around  zero  by  lag  B  or  so.  This  gives  us  298x16,  or  4768  numbers 
which  should  behave  as  an  independent  identically  distributed  sample  of  Gaus¬ 
sian  random  variables.  These  numbers  had  a  kurtosis  of  3.167  and  a  skew  of 
.058  respectively  A  test  against  the  hypothesis  of  kurtosis  exactly  3  would 
reject  at  the  2%  level;  however  a  similar  test  for  zero  skew  is  only  significant  at 
about  the  10%  level.  Considering  how  sensitive  sample  kurtosis  is  to  tails  larger 
than  expected  from  normality,  this  result  is  remarkably  close  to  the  normal 
theory  prediction. 

QQ  plots  provide  a  very  good  qualitative  measure  of  Gaussianity-  deviations 


from  the  expected  straight  line  are  easy  to  see.  This  plot  is  shown  in  Exhibit  4.2. 
Qualitatively,  the  fit  to  a  straight  line  seems  excellent.  Quantitatively,  we  can 
test  this  fit.  This  test  is  equivalent  to  a  Kolmogorov-Smirnov  test  of  whether  the 
distribution  is  Gaussian  and  only  rejects  the  hypothesis  that  the  distribution  is 
Gaussian  at  a  significance  level  of  greater  than  20%  Thus  the  subsample  of  4768 
points  passes  this  test  also.  A  histogram  of  the  subsample  with  a  superimposed 
exact  normal  curve  is  shown  in  Exhibit  4.3.  This  picture  also  qualitatively 
confirms  the  hypothesis  of  Gaussianness. 

We  conclude  that  the  assumption  of  a  Gaussian  distribution  for  this  stretch 
of  the  data  appears  well  substantiated,  except  for  three  records  containing 
bursts.  It  thus  seems  a  reasonable  hypothesis  that  the  background  noise 
remains  Gaussian,  even  in  the  high  variance  regions,  especially  if  we  can  demon¬ 
strate  that  the  background  noise  keeps  much  the  same  power.  This  is  one  of  the 
conclusions  in  the  next  section. 

4.1.2.  Existence  and  Identification  of  Tones. 

In  the  trace  of  record  variances,  there  are  several  high  variance  regions: 
300-1400,  1600-1900  and  3300-4300  (regions  2,  3,  5  and  7).  Minor  fluctuations 
are  apparent  in  regions  2600-2800,  and  perhaps  5000-5400. 

A  reasonably  typical  log  periodogram  (smoothed)  of  a  group  of  these 
records  (3508-3512)  is  pictured  in  Exhibit  4.4.  Notice  there  is  an  obvious  double 
spike  at  about  1250  Hz  and  a  possible  harmonic  of  these  spikes  at  2500  Hz.  So 
tones  occur  some  of  the  time.  Inspection  of  periodograms  in  the  first  group 
(300-1300)  reveal  a  more  complicated  structure  (e  g.  see  Exhibit  4.5,  record 
730).  but  we  suspect  that  the  high  variance  records  observed  in  this  group  have 
the  same  cause  We  shall  treat  this  group  separately,  later  in  this  section. 

We  first  note  one  possible  consequence  of  Model  (1)  The  addition  of  a 
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discrete  tone  adds  a  sinusoidal  component  to  the  observed  data.  If  the  tone  has 
large  power  compared  to  the  background  noise,  then  the  distribution  of 
observed  series  will  be  similar  to  the  distribution  obtained  by  randomly  sam¬ 
pling  X  from  a  distribution  with  distribution  function  G,  where 

G(x)  =  ~r+  ^-Arcsin  x 
'  2  rr 

The  kurtosis  of  this  distribution  is  1.5.  It  is  to  be  expected  that  the  addition  of  a 
small  number  of  tones  with  high  power  will  have  much  the  same  effect,  namely 
lowering  the  observed  kurtosis  to  less  than  3.  This  need  not  necessarily  be  true; 
if  there  are  harmonics  with  power  close  to  that  of  the  fundamental  tone,  the 
kurtosis  may  vary  over  a  range,  including  values  larger  than  3. 

This  suggests  that  a  comparison  of  low  kurtosis  records  with  their 
corresponding  variances  might  be  revealing.  These  plots  (Exhibit  3.8  and  Exhi¬ 
bit  3.2)  clearly  pick  out  the  regions  mentioned  above,  except  for  5300-5500  and 
perhaps  2600-2800.  The  regions  1600-1900  and  3300-4300  are  particularly  strik¬ 
ing,  with  almost  no  kurtoses  greater  than  three.  The  case  of  few  tones,  with  few 
harmonics  seems  to  be  true  here.  However,  the  region  300-1300  retains  many 
records  with  kurtoses  greater  than  3. 

Note  that  this  does  not  prove  Model  (l),  even  in  the  regions  of  low  kurtosis. 
Strengthening  the  case  for  Model  (i)  requires  a  demonstration  that  the  back¬ 
ground  noise  is  not  varying,  and  that  larger  kurtoses  in  the  300-1300  record 
number  region  are  actually  due  to  tones  with  many  harmonics,  and  not  to  other 
reasons. 

It  is  certainly  not  unreasonable  to  expect  the  background  noise  to  change 
at  the  same  time  as  we  observe  mixtures  of  discrete  frequencies  (a  failure  of 
Model  (1)).  But,  as  discussed  in  the  beginning  of  Section  4,  the  background 
noise  appears  very  constant  with  virtually  no  change  in  power  for  long  record 


stretches. 


-  16  - 


To  check  this  we  wished  to  find  the  power  due  to  the  tones  and  remove  it, 
leaving  only  the  power  due  to  the  background  noise.  We  did  this  by  smoothing 
the  periodogram  in  such  a  way  as  to  remove  the  influence  of  tones  at  or  around 
discrete  frequencies  while  preserving  the  overall  shape  of  the  periodogram 
Since  the  total  variance  or  power  of  a  process  is  just  the  area  under  the  periodo¬ 
gram,  we  could  then  compute  the  "residual  variance"  of  the  part  of  the  noise 
modelled  by  h  by  taking  the  area  under  the  smoothed  periodogram.  Such  a 
smoothing  procedure  should  have  the  property  that  it  preserves  the  shape  of 
periodograms  from  records  with  no  apparent  spikes  due  to  tones.  Classical 
smoothing  methods  Eire  not  applicable,  as  they  merely  average  the  spikes  at 
discrete  frequencies  into  the  neighboring  frequencies  In  our  circumstance, 
where  the  variance  attributable  to  a  tone  is  between  two  or  three  times  the  vari¬ 
ance  at  all  other  frequencies,  classical  smoothing  will  clearly  produce  a  highly 
biased  estimate  of  the  power  due  to  h* . 

The  usual  method  of  removing  tones  is  to  use  complex  demodulation  to 
identify  the  frequency  at  which  the  tone  occurs,  and  remove  that  frequency  with 
some  sort  of  average  of  surrounding  frequencies.  The  major  problem  with  this 
approach  is  that  it  does  not  lend  itself  well  to  automation,  for  it  requires 
automatic  identification  of  the  frequencies  at  which  tones  occur.  Alternatively, 
one  can  attempt  to  use  some  sort  of  sequential  testing  procedure  like  the  P(\) 
test  Priestley  (1981)  advocates,  and  then  remove  the  peaks  A  double  windowing 
teahnique  like  that  suggested  by  Hannan  (1961),  is  also  a  possibility  All  these 
approaches  require  considerable  work  which  is  not  easy  to  automate. 

Some  sort  of  non-linear  smoother  seems  to  answer  these  problems  A  non¬ 
linear  smoother  will  produce  a  biased  estimate  of  h,  but  under  the  assumption 
of  sufficiently  smooth  h,  and  Gaussianity  in  the  background  noise,  this  bias  can 
be  adjusted  for  reasonably  well.  The  smoother  we  chose  to  use  was  that  of  a 
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(2m  +  1)  point  running  median.  Denote  this  smoother  by  Sm  Then  the  value  of 
the  smoothed  periodogram  at  Fourier  frequency  c*  =  2m/ N  is: 

SmP(Ui )  =  median  $/,(uj){ 

We  need  an  endpoint  rule:  we  chose  to  simply  left  extend  P  for  values  of  i<m; 
similarly  for  right  extension.  One  of  the  virtues  of  this  smoother  is  that  it  gives 
the  same  smooth  regardless  of  whether  one  operates  on  the  raw  periodogram 
and  takes  logs,  or  operates  with  the  log  periodogram,  since  log  is  a  monotone 
transformation. 

Intuitively,  this  smoother  will  remove  any  narrow  spikes  in  the  periodo¬ 
gram.  As  an  estimator  of  h  it  has  several  sources  of  bias.  Firstly,  it  will  smooth 
out  peaks  and  troughs  in  h.  That  this  bias  is  kept  small  will  depend  on  A  being 
sufficiently  smooth  in  relation  to  the  ’’peakiness"  of  the  spikes  in  the  periodo¬ 
gram  due  to  tones.  Secondly,  even  if  h  is  quite  smooth,  SmP  will  have  a  con¬ 
sistent  bias  due  to  the  fact  that  the  median  of  the  exponential  distribution  is 
log(2)  rather  than  1.  This  can  easily  be  adjusted  for.  Thirdly,  it  will  fail  if  there 
are  too  many  spikes  too  close  together  in  the  periodogram.  In  this  case  it  is 
unclear  how  to  mechanize  a  procedure  for  tone  removal  —  this  appears  to  be  a 
problem  in  the  300-1300  group,  as  discussed  a  little  later. 

The  major  theoretical  justification  for  using  Sm  can  be  found  in  Theorem  1 
of  Appendix  A,  which  states  that,  under  mild  conditions  on  h  and  m, 
Sn(Pf,(d))/ \og(?)  converges  to  h(d)  uniformly  in  probability.  Further,  uniform 
convergence  to  h  in  probability  implies  that  the  residual  variance  given  by 

4rS  sm(P(u))  -*  fh(u)du  =  var(A')  in  prob 

yv  i=i 

where  X  is  the  Gaussian  part  of  the  time  series.  (For  am  exact  statement,  see 
Corollary  1  in  Appendix  A). 

We  chose  m=25  This  removes  the  discrete  spikes  almost  totally,  and,  since 
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h  looks  locally  smooth,  the  residual  variance  should  not  have  too  much  bias  in 
it.  Model  (1)  predicts  that  the  residual  variance  trace  should  be  constant.  We 
calculated  this  residual  variance  for  each  record  and  the  result  of  our  pro¬ 
cedure  can  be  seen  in  Exhibit  4.6,  where  the  original  sample  variance  trace  is 
plotted  above  the  residual  variance  trace.  (All  records  containing  suspected 
bursts,  i.e.  those  with  sample  kurtoses  bigger  than  4,  are  not  plotted  —  our 
smoothing  procedure  can  only  handle  tones,  not  bursts.)  The  procedure  seems 
to  do  a  reasonable  job  of  normalizing  the  variance  trace  everywhere  except  for 
the  region  of  records  300-1000.  That  these  records  still  fit  into  the  framework  of 
Model  (1)  remains  to  be  shown.  It  is  clear  though,  that  Model  (1)  is  only  approxi¬ 
mately  true  elsewhere. 

The  region  300-1000  is  anomalous  in  other  ways.  The  kurtoses  are  not  con¬ 
sistently  low.  An  inspection  of  the  skew  trace,  with  skew  values  limited  to 
between  +.6  and  -.6  (Exhibit  3.8)  clearly  indicates  that  this  region  is  qualita¬ 
tively  different  from  the  region  3300-4300,  where  skew  appears  to  decrease  in 
absolute  value.  Record  730  is  displayed  in  Exhibit  4  5,  and  clearly  shows  a 
strong  non-sinusoidal  but  periodic  component.  The  periodogram  of  this  record 
shows  the  fundamental  at  about  240  Hz  and  six  or  seven  harmonics,  five  of  which 
have  about  the  same  power  The  running  median  smoother  cannot  deal  with  this 
case,  as  the  Fejer  kernels  about  each  peak  leak  substantial  amounts  of  power 
which  reinforce  each  other,  at  least  for  the  length  of  running  median  we  chose. 
It  is  obvious  that  for  a  sample  length  of  1024,  no  running  median  smoother  will 
work  well  in  this  case.  Other  periodograms  we  looked  at  have  a  yet  more  com¬ 
plex  structure  with  two  or  three  fundamentals  apparent,  together  with  their  har¬ 
monics 

We  believe  that  Model  (1)  is  still  valid  for  the  whole  region,  but  this  is  on 
purely  observational  grounds  -  the  periodograms  we  considered  look  like  they 


retain  the  same  underlying  structure.  This  remains  an  area  where  we  need  to 
do  more  rigorous  research. 

There  remains  open  the  question  of  how  one  might  modei  the  frequencies 
and  amplitudes  of  these  tones,  as  a  random  process.  This  is  not  an  easy  task,  as 
it  seems  apparent  that  the  tones  at  different  frequencies  are  connected  with 
each  other  -  to  wit,  it  is  clear  that  we  sometimes  observe  a  fundamental  with  its 
harmonics.  Any  reasonable  model  must  take  such  structure  into  account.  Con¬ 
sideration  of  the  physics  ef  underwater  acoustics  and  of  acoustic  energy 
released  in  the  fracture  of  solids  might  suggest  some  avenues  of  approach. 

That  these  tones  have  a  complex  structure  even  without  considering  the 
dependence  structure  at  one  time  is  illustrated  by  the  fact  that  even  in  periods 
of  a  second  or  so,  the  tone  may  wander  in  both  frequency  and  in  amplitude.  A 
plot  of  the  spikes  in  the  estimated  spectra  of  a  series  of  records  1770-1830  (see 
Exhibit  4.7)  illustrates  this  point.  This  plot  was  made  by  first  computing  a  spec¬ 
trum  estimate  for  each  of  the  records.  The  estimate  used  was  a  length  five 
Daniell  window  applied  to  the  periodogram.  Then  the  difference  between  this 
estimate  and  a  smoothed  version  (smoothed  using  running  medians  of  length  41) 
was  computed  in  order  to  leave  only  the  spikes  The  plot  shows,  for  each  record, 
only  the  upper  15%  of  these  differences. 

In  the  early  records,  1770  to  1790  or  so,  there  is  quite  a  strong  tone  at 
around  1100  or  1200  Hz,  which  appears  to  split  into  two  tones;  one  staying  at  the 
same  frequency  until  late  in  the  sequence,  the  other  wandering  up  to  about  1500 
Hz  The  two  tones  seem  to  be  merging  again  by  record  1830.  There  is  also  an 
apparent  harmonic  at  2200-2300  Hz,  which  increases  in  frequency  just  as  its  fun¬ 


damental  does. 
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5.  Summary 

This  report  has  modeled  the  undersea  noise  with  a  model  suggested  by  an 
exploration  of  the  dataset  itself.  The  model  consists  of  three  parts;  the  first 
part  is  a  relatively  stable  Gaussian  process  (called  the  background  noise  in  this 
report)  which  is  assumed  to  have  a  smooth  continuous  spectrum.  Superimposed 
on  this  Gaussian  process  are  a  small  number  of  random  tones.  The  frequencies 
amplitudes  of  these  tones  are  changing  over  time  and  are  sometimes  absent 
entirely.  The  third  component  in  the  model  is  the  appearance  of  sporadic, 
extremely  transient,  high  intensity  bursts  of  impulsive  sound  energy  (bursts) 
These  bursts  are  quite  rare  in  the  data  we  examined. 

The  model  we  suggest  is  incomplete  in  several  respects.  It  is  only  approxi¬ 
mate,  but  we  believe  this  is  not  a  problem.  This  belief  could  best  be  tested  by 
trying  out  the  model  with  another  dataset.  The  model  is  unproven  in  some 
regions  of  the  data,  particularly  where  a  large  number  of  tones  are  difficult  to 
distinguish  from  the  background  Gaussian  noise. 

There  are  several  areas  remaining  open  for  investigation.  One  such  area  is 
in  further  modeling.  We  do  not  specifically  estimate  the  parameters  associated 
with  the  Gaussian  process.  In  addition,  we  do  not  attempt  to  model  any  proba¬ 
bilistic  process  for  the  impulsive  noise  bursts,  nor  do  we  attempt  to  modei  the 
tones  by  a  probability  model.  Trying  to  model  the  bursts,  perhaps  as  Poisson 
events,  may  be  fruitful.  The  problems  associated  with  modeling  the  tones 
appear  more  difficult.  For  one  thing,  we  cannot  always  dissect  out  the  tones,  as 
mentioned  earlier,  and  for  another,  the  tones  probably  have  a  complex  struc¬ 
ture  dependent  on  the  physics  of  acoustic  energy  released  due  cracking  and 
strains  in  the  ice  However,  it  may  be  that  a  simple  model  useful  in  signal  detec¬ 
tion  can  adequately  describe  these  tones. 

Another  area  of  interest  is  in  signal  detection  Knowledge  of  the  noise  can 
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suggest  useful  methods  of  signal  detection.  There  are  a  variety  of  approaches  -- 
for  example,  one  might  derive  theoretical  results  which  can  be  tested  by  adding 
signal  to  the  noise  data.  This  might  best  be  done  in  conjunction  with  more  com¬ 
plete  modeling  Another  interesting  approach  would  be  to  compare  the 
effectiveness  of  existing  detection  schemes  in  the  presence  of  this  noise. 


Appendix  A 


The  aim  of  this  appendix  is  to  give  the  technical  results  used  in  the  paper. 
These  are  given  in  Theorem  2  and  its  Corollary. 

Suppose  that  we  observe 

Zi=Xi  +  Yi.i  =  1.2 . N,  (1) 

where  is  Gaussian  with  non-normalized  spectral  density  h{d),  positive  for  all 

cj,  and  Yx  is  a  mixture  of  sinusoids  at  I  discrete  frequencies  with  amplitudes 
and  arbitrary  phase  shifts  (i  =  1,  2 .  /).  Assume  that  Xx  is  a  linear  pro¬ 

cess: 

Xi=tPtti-t  (2) 

where  the  are  independent  and  Gaussian,  E(st)  -  0.  E{tf)  =  1,  and  the  ft  are 
constants  satisfying  the  condition 

ijc  (3) 

for  some  c  >  0.  Condition  (3)  implies  that  the  spectral  density  h  is  continuous 

and  bounded  on  [O.ttJ.  Lastly,  assume  h{6)  >  0,  for  ail  u. 

Form  the  standard  periodogram  of  Zx  at  Fourier  frequencies  Uj n  =  Znj/ N , 
0  <  j  <  N/  2  (we  ignore  the  end  points  0,  N /  2  so  that  Theorem  1  holds  true).  So 

pN(ujN)  =  (o 

We  need  two  preliminary  results-  the  first  is  a  lemma  due  to  Hoeflfding 
(1963): 

Lemma  1 

Let  Wx,  i  =  l,  2 m  be  independent  random  variables  bounded  between  zero 

and  one  Then,  for  each  t  >  0, 
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The  second  result  is  due  to  Walker  (1965)  and  concerns  the  asymptotic 
independence  of  Pf/{u}f/)  for  different  Fourier  frequencies: 

Theorem  1 

Let  Xi  be  given  by  (2)  and  satisfy  (3).  Then: 

PnWjn)  ~  IjjV  +/?A’(ityv)  (5) 

where  the  ViN  are  independent,  each  with  distribution  4-h(^jjv)  x!  (Notice  the 

Fourier  frequencies  do  not  include  0,  N /  2). 

Further: 

0<m^2  =  OpC)  (6) 

Under  these  conditions  we  prove  that  the  running  median  estimate  of  h 
converges  uniformly  in  probability  by  proving  the  following  theorem. 

Theorem  2 

Let  Zt  be  a  process  as  given  in  (l).  satisfying  condition  (3).  Let  m  be  an  integer 
depending  on  TV  such  that  as  TV-**,  the  following  holds: 

-*  0  as  TV  -*  *  (7) 

and  there  exists  a  y  >  ^-so  that: 

—p  -*  °°  as  TV  -»  *  (6) 

and  assume  h  bounded  away  from  zero. 

Given  u,  let  [i)f  {  be  the  set  of  m  oiN' s  nearest  to  u  (mod  n).  Let 

K^n{»)  =  median  of  |Pv(i>;-)/  log(2)j  (9) 

Let  t  >  0  be  given.  Then,  as  TV-**, 

P[supihmAr(cj)-/i(cj)':  >  e]  -*  0  (10) 

Proof 

First  note  that  Pn{U)n)  are  approximately  distributed  as  TV  independent 
(ie.  exponentials  with  parameter  but  with  a  fixed 
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constant  added,  g(u jx),  where 

0(a)  =  y£lA?\lrN(u~Si)+FN(u+?l)  j  (II) 

*<= 1 

where  the  /Vs  and  <Pi ’s  come  from  the  discrete  frequencies  assumed  in  the 
theorem  and  FN  is  a  Fejer  kernel; 


Fn<$)  =  _j_  sin2(Ayi?/  2) 
2  nN  sin8(tf/2) 


(12) 


Given  the  decomposition  in  Theorem  1,  we  know  the  distribution  functions  of  Vjjf, 
say  they  are  given  by  Gjn  For  simplicity  in  the  following  argument,  let  us 
further  suppose  7=1.  The  case  of  7>1  is  proved  along  the  same  lines,  except  for 
increased  notation.  Fix  an  u0.  Let  the  median  of  an  exponential,  parameter 
7i(uc)-1,  be  given  by: 


io  =  log(2)./i(«c) 

Let  £m,v'(cJc)  =  log(2).hmjv(iJo)  We  shall  abuse  notation  by  dropping  the  argument 
to  Cton'  when  it  is  clear  to  which  we  are  referring.  Then  tm#'  is  an  estimate  of 
£0.  the  median  of  an  exponential,  parameter  A(ac).  We  shall  work  with  £0  and  the 
estimate  so  without  loss  of  generality,  replace  e  by  e.iog(2)  in  the  state¬ 
ment  of  this  theorem. 

Make  the  decomposition  (5).  We  work  first  with  the  case  of  considering 
£»n/v(uc)  =  the  median  of  j  Vfn\  and  extend  to  later.  Again,  we  drop  the  argu¬ 
ment  to  £rnx  if  it  is  clear  which  argument  is  implicitly  referred  to.  So: 

F l  I  £mA'  — £0  >  e]  =  P[^mN  ^  fc  "*■  P[$mN  ^  fo  —  E]  (»3) 

Let 


Then 


a  =  f0  +  f 


P[tmN  >  a]  =  P[  g-or  more  PN{uf)  >  a] 

=  P[2/[Pw(Q;-)>a]^  f-] 

t*i  c 

=  P[£V[P*W)  >  a]-GfM)*f,{±—GfN(a))} 

<■1  i=i  c 


=  exp(-2m^v) 


by  Lemma  1,  where 

<5mA'  =  jj-) 

We  now  show  for  m,  N  large  enough 

<5Lv  >  0 

Put  L  =  maxf/i(iO_I;i>  in  [0,27r)j.  L  exist  by  the  assumption  h  is  positive.  As 
maximum  of  h  is  bounded,  put  M  =  min{/i(y)_1J,  and  for  N  large  enough, 

A\/  N  <  for  all  u 

So  assume  this  holds.  Also: 

AiFf/iojn-?')  s  h(cjjAr)-1  for  at  most  k\  of  the  ^ 

N 

We  can  choose  fc.v  to  satisfy  -~r-  -*  0.  where  fc*  comes  from  consideration  of  the 

v 

form  of  the  Fejer  kernel,  which  is  bounded  by  for  the  Fourier 

points  farthest  from  <pi. 

Simultaneously,  we  can  choose  fcy  to  satisfy - »  0.  by  considering  condition 

m 

(8);  e  g  choose  *»/  =  N*.  where  ^-<  {  <  y. 

So.  for  at  most  k #  points,  fo+e-3  (“/)  <  0  So: 

6mit  =  G^(a)_r 

m  j? 

except  for  at  most  kN  cj;-, 

=  2^£r.i-exp{-^(Wj  )",c+[ft(wc)-1-^("/)'1]fo+^(«j')"13(«j )!]  (15) 

except  for  at  most  kN  uf .  We  use  now  uniform  continuity  of  h-1  and  the  fact 
m/  N  -*<*>  to  note  that  for  m  large  enough, 


IMu')"1  ~h(»)-'\  <e'<7T 

4so 


if  I  cjf  — cj0  !  <m/N.  So  (15)^ 


*  ikf  I1  'exp(~l£  +  ^ +  "*/JV>]  -  tr 

where  we  have  dealt  with  the  "at  most  kN  Uj-  's"  by  bounding  any  exceptional 
term  G}-s  (a)-4-from  below  by  -1/  2. 


*1  -Sr>«o>o 


for  m,  N  large  enough,  recalling - *  0.  It  is  clear  that  this  bound  does  not 

TTL 

depend  on  cj0  So  we  have  proved 

P[tmN  >  lo  +  e]  *  b  'Zmif 

for  m,  A'  large  enough.  We  prove  an  exactly  similar  bound  on  the  second  term 
in  the  right  side  of  (13).  Let  6X  be  the  minimum  of  the  two  50's.  Hence 

P[\Uv~to\  >e]*2a~*m>f  (16) 

Now  we  bound  the  effect  of  including  PyiVjx): 

P  [SUD  !£nVjv,(u0Wol  >  £J  =  PL  !  >  c] 

Ofiu^tn  0<j<N/2  J  J 

(where  =  log(2 )  h{ujN)) 

*  I  lN(UjN)-*j  !  >  E  -SUp  |  RN{Uj  at)  |  |  J 

*  P:0<®2!^(^-v)-^  1  >  S/21 

+  ^[0<SU^/z!^(^;w)i  >  e/2]  (17) 

The  second  term  of  the  (17)  is  o,(l)  by  (6).  We  show  that  the  first  term  is  also 
o,(l).  But  by  applying  (16), 

pL0<i<l/2!frnJv(^vv)_fJ !  >£/2]-  7Ve_2m{Z 
(by  applying  Lemma  1  and  sub-additivity  of  probability  measures) 

= exp( -2m  <5 2 + log  AO 
„  .,  m 

-*  0  if - >  00 


This  is  implied  by  condition  (5),  so  we  are  done. 

We  have  an  immediate  corollary  to  Theorem  1: 

Corollary  1 

Let 

J'mA'  =  ~rr S 

/v  i  =  l 

and  let 

V  =  fh{u)du  =  var(A") 

Then,  given  e  >  0,  and  under  the  assumptions  of  Theorem  A. 

P[\ymN-V\  >e]-0 

as  N  -*  ». 

Proof 

The  proof  is  immediate  on  considering  the  result  of  the  theorem  and  recalling 
that  h  is  uniformly  continuous. 
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